System and method for high-order accurate device model approximation

ABSTRACT

An system and method are disclosed for efficiently approximating analytical circuit device models. A preferred embodiment includes a method for obtaining smooth and accurate approximations of analytical device models, comprising the steps of identifying a first set of measurement units; locating two or more sets of units that neighbor one or more of said measurement units; for each set of the two or more sets of neighbor units, obtaining the union of one or more of said sets of neighbor units and the first set of measurement units; calculating the smoothness of the analytical device model within one or more of said unions; and selecting at least one of said unions within which the analytical device model is the smoothest as the new set of measurement units.

FIELD OF THE INVENTION

[0001] The present invention relates generally to the field of devicemodel approximation. More particularly, the present invention relates toan efficient system and method for approximating analytical circuitdevice models.

BACKGROUND OF THE INVENTION

[0002] Simulation is an indispensable verification step beforecommitting an integrated circuit to an expensive manufacturing process.An important step in circuit simulation is model evaluation. Modernanalytical models such as bsim3 and bsim4 have become more and morecomplicated and expensive to evaluate in simulation. The percentage ofmodel evaluation time in circuit simulation can be as high as 70-80%.This number can grow to 90% if one includes simulation time spent inusing all the individual device evaluations to generate the circuitequations. As a result, model evaluation has become a bottle neck in theimprovement of circuit simulation efficiency.

[0003] It is becoming more and more difficult to maintain some degree ofsmoothness in the approximated analytical model due to its complexity.As a result, many resources have been spent in checking the smoothnessof device models and the consistency between function values and theirderivatives in both the Electronic Design Automation (EDA) and thesemi-conductor industries.

[0004] A high-order accurate table model is a potential solution to theabove problems. Table models provide polynomial approximations ofanalytical models based on table-grid interpolations of the originalcurrent-voltage and charge-voltage (i-ν/q-ν) curves. Most existing tablemodels use spline interpolation based on a fixed stencil consisting of afixed group of consecutive interpolation points. Fix stencilinterpolations work well if the curve to be approximated is globallysmooth. The problem associated with fixed stencil interpolations,however, is that for curves containing a C⁰ corner point (i.e. curvesthat are continuous but not differentiable at the comer point), fixedstencil interpolation of second or higher order accuracy will benecessarily oscillatory near the C⁰ comer point.

[0005] The impact of this oscillatory behavior can be illustrated usingthe example in FIG. 1. Shown in FIG. 1 is the 3^(rd)-order polynomialapproximation with a fixed stencil for the step function. Notice theobvious undershoots near the bottom corner point. Such oscillations aredenoted as Gibbs Phenomena in spectral methods. The main characteristicsof these oscillations is that they do not decay in magnitude when themesh or table grid is refined. In other words, no matter how big thetable size is (or how large the number of interpolations is), theoscillations still remain at the same magnitude. Besides causinginaccuracy, the non-diminishing oscillations may cause the Newtoniteration method to fail in solving nonlinear circuit equations.

[0006] Besides causing non-diminishing oscillations, existing tablemodel-based simulation methods also lead to unnecessary memory usage andinefficient memory access. This is because these methods construct andstore all the measurement data entries (i. e. table entries) at thebeginning of simulation. Storing all the data entries is not necessarysince simulations usually only concentrates on “local” data entriessurrounding one or more operating points. Moreover, these data entriesmay contain redundant information since boundaries of these data entriesmay overlap. Storing these redundant information is certainlyunnecessary. Moreover, storing all the table entries at the beginning ofsimulation makes it impractical to build the table entries adaptivelyand accurately due to the huge memory consumption.

[0007] There is therefore a need to develop a method and system forobtaining smooth, accurate, and computationally efficient approximationof analytical device models.

SUMMARY OF THE INVENTION

[0008] The presentation invention includes a method for obtaining smoothand accurate approximations of analytical device models, comprising thesteps of identifying a first set of measurement units; locating two ormore sets of units that neighbor one or more of said measurement units;for each set of the two or more sets of neighbor units, obtaining theunion of one or more of said sets of neighbor units and the first set ofmeasurement units; calculating the smoothness of the analytical devicemodel within one or more of said unions; and selecting at least one ofsaid unions within which the analytical device model is the smoothest asthe new set of measurement units.

[0009] The present invention also includes a method for obtaining smoothand accurate approximations of analytical device models in a(m+1)-dimensional manifold where m is an integer, comprising the stepsof identifying a first set of measurement units in the (m+1)-dimensionalmanifold; locating two or more sets of units that neighbor one or moreof said measurement units; obtaining, for each m-dimensional boundary ofthe (m+1)-dimensional manifold, an approximation of the analyticaldevice model based on the two or more sets of neighboring points and thefirst set of measurement units; and blending the approximations obtainedfor each m-dimensional boundary of the(m+1)-dimensional manifold toproduce an approximation of the analytical device model.

[0010] The present invention further includes a method for reducingmemory usage during circuit simulations, comprising the steps ofidentifying one or more operating points; obtaining, for each of the oneor more operating points, a set of measurement points surrounding saidoperating point; eliminating overlapping boundaries among each set ofmeasurement points; and constructing a hash table for storing each setof measurement points.

BRIEF DESCRIPTION OF THE DRAWINGS

[0011]FIG. 1 compares the approximations of a step function obtainedusing the essentially non-oscillatory (ENO) approximation method and aprior art approximation method.

[0012]FIG. 2 schematically depicts a flow chart of the ENO-basedapproximation method.

[0013]FIG. 3 schematically depicts a flow chart of an embodiment of theENO-based approximation method in a m-dimensional manifold.

[0014]FIG. 4 schematically depicts the table data storage scheme.

[0015]FIG. 5 schematically depicts the non-uniform grid table.

[0016]FIGS. 6-10 demonstrate simulation results using the ENO-basedapproximation method.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0017] One preferred embodiment of the present invention uses dynamicstencils to obtain smooth and accurate approximations of analyticaldevice models. A flow chart of this embodiment is depicted in FIG. 2. Atstep 210, a first set of measurement units is identified. Step 220locates two or more sets of neighboring units that neighbor one or moreof said measurement units. At steps 230-250, for each set of the two ormore sets of neighboring units, the dynamic stencil scheme obtains theunion of said each set and said first set of measurement units andcalculates the smoothness of the analytical device model within theunion. At step 260, the dynamic stencil scheme selects at least one ofsaid unions within which the analytical device model is the smoothest asa new set of measurement units. Steps 220-260 are repeatedly iterateduntil the new set of units achieves the desired accuracy of theanalytical device model. The approximated analytical device model isobtained at step 280 based on the final set of measurement units and thesmoothness calculations.

[0018] The above approximation process can be performed in any subspaceof the measurement space and the corresponding measurement units in thesubspace can be of any shape. In one embodiment, the approximationprocess further comprises, before step 210, the steps of providing apartition of the measurement space, and for each partition, executingsteps 210-260. The approximated analytical device model is obtained atstep 280 based on the final set of measurement units for each partition.

[0019] Examples of the dynamic stencil scheme include the EssentiallyNon-Oscillatory (ENO) schemes, which are described by Harten et al. in[1], the contents of which is hereby incorporated in its entirety byreference. In one embodiment, the original charge (q) curve describingthe analytical device model is described by the mapping q: V→R, where Vis a m-dimensional manifold representing the measurement index and R isthe set of real numbers representing the measurement value of theoriginal curve. Measurement units and neighboring units are identifiedbased on a measurement index contained in manifold V. The measurementindex can be of any shape. In one embodiment, this index is a gridconsisting of cells of equal or variable size. Measurement units andneighboring units thus obtain their value from selected grid-points onthe grid.

[0020] The dynamic stencil grows based on the smoothness calculation ofthe original curve within sets of stencil points. In one embodiment, thesmoothness of the original curve may be obtained by calculating theNewton divided difference of the original curve within the sets ofstencil point. Specifically, when the index grid is one-dimensional, themeasurement units become the “stencil points”. The smoothnesscalculation and the selection of stencil points can be described asfollows:

[0021] Given a grid${a = {v_{\frac{1}{2}} < v_{\frac{3}{2}} < \quad {\ldots \quad v_{N - \frac{1}{2}}} < v_{N + \frac{1}{2}}}},$

[0022] define cells, cell centers, and cell sizes by${I_{i} \equiv \left\lbrack {_{i - \frac{1}{2}},_{i + \frac{1}{2}}} \right\rbrack},{_{i} \equiv {\frac{1}{2}\left( {_{i - \frac{1}{2}} + _{i + \frac{1}{2}}} \right)}},{{\Delta \quad _{i}} \equiv {_{i + \frac{1}{2}} - _{i - \frac{1}{2}}}},{i = 1},2,\ldots \quad,{N.}$

[0023] Denote Δν as the maximum cell size and represented it byΔν≡max_(I≦i≦N)Δν. The approximation problem can be defined as:

[0024] Given the cell average of the original curve q(ν):${{\overset{\_}{q}}_{i} \equiv {\frac{1}{\Delta \quad _{i}}{\int_{_{i - \frac{1}{2}}}^{_{i + \frac{1}{2}}}{{q(\xi)}\quad }}}},{i = 1},2,\cdots \quad,N,$

[0025] find a polynomial p_(i)(ν), of degree at most k−1, for each cellI_(i), such that it is a k-th order accurate approximation to thefunction q(ν)inside I_(i), i.e.

p _(i)(ν)=q(ν)+O(Δν^(k)), νεI _(i) ,i=1, . . . , N.  (1)

[0026] This gives approximations to the function at cell boundarieswhich are also -th order accurate, i.e. $\begin{matrix}{{{p_{i}\left( _{i + \frac{1}{2}} \right)} = {{q\left( _{1 + \frac{1}{2}} \right)} + {O\left( {\Delta \quad v^{k}} \right)}}},{{p_{i}\left( v_{i - \frac{1}{2}} \right)} = {{q\left( v_{i - \frac{1}{2}} \right)} + {O\left( {\Delta \quad v^{k}} \right)}}},{v\quad ɛ\quad I_{i}},{i = 1},\cdots \quad,{N.}} & (2)\end{matrix}$

[0027] Function p_(i)(ν)can be replaced by other simple functions, suchas exponential functions [2]. For large-scale circuit simulation,polynomial may be the most efficient to evaluate.

[0028] The 0-th degree Newton divided differences of function aredefined by${{Q\left\lbrack v_{i - \frac{1}{2}} \right\rbrack} \equiv {q\left( v_{i - \frac{1}{2}} \right)}},{i = 1},2,{\cdots \quad {N.}}$

[0029] In general, the j-th degree Newton divided differences for j≧1,are defined inductively by${Q\left\lbrack {v_{i - \frac{1}{2}},\cdots \quad,v_{i + j - \frac{1}{2}}} \right\rbrack} \equiv \frac{{Q\left\lbrack {v_{i + \frac{1}{2}},\cdots \quad,v_{i + j - \frac{1}{2}}} \right\rbrack} - {Q\left\lbrack {v_{i - \frac{1}{2}},\cdots \quad,v_{i + j - \frac{3}{2}}} \right\rbrack}}{v_{i + j - \frac{1}{2}} - v_{i - \frac{1}{2}}}$

[0030] The Newton form of the k-th degree interpolation polynomial P(ν),which interpolates q(ν) at k+1 points starting from grid point$v_{i - r - \frac{1}{2}}$

[0031] can be expressed using the divided differences by${P(v)} = {\sum\limits_{j = 0}^{k}\quad {{Q\left\lbrack {v_{i - r - \frac{1}{2}},\cdots \quad,v_{i - r + j - \frac{1}{2}}} \right\rbrack}{\prod\limits_{m = 0}^{j - 1}\quad \left( {v - v_{i - r + m - \frac{1}{2}}} \right)}}}$

[0032] We can take the derivative of P(ν)to obtain p(ν):${p(v)} = {\sum\limits_{j = 0}^{k}\quad {{Q\left\lbrack {v_{i - r - \frac{1}{2}},\cdots \quad,v_{i - r + j - \frac{1}{2}}} \right\rbrack}{\sum\limits_{{m = 0},{m \neq r}}^{j - 1}\quad {\prod\limits_{m = 0}^{j - 1}\quad \left( {v - v_{i - r + m - \frac{1}{2}}} \right)}}}}$

[0033] An important property of divided difference is:$\left. {{Q\left\lbrack v_{i - \frac{1}{2}} \right\rbrack},\cdots \quad,v_{i + j - \frac{1}{2}}} \right\rbrack = \frac{q^{(j)}(\xi)}{j!}$

[0034] for some ξ inside the stencil$v_{i - \frac{1}{2}},{< \xi < v_{{i + j - \frac{1}{2}},}},$

[0035] as long as the function q(ν) is smooth inside the stencil. Thusthe divided difference is a measurement of the smoothness of the curveinside the stencil.

[0036] Based on these definitions, the dynamic selection of stencils canbe described as follows. Suppose we want to find a stencil of k+1consecutive points, which must include${v_{i - \frac{1}{2}}\quad {and}\quad v_{i + \frac{1}{2}}},$

[0037] such that q(ν)is the “smoothest” in this stencil comparing withother possible stencils. This job may be performed in the preferredembodiment by the following steps. In each step only one point may beadded to the stencil. Start with the two point stencil${{\overset{\_}{S_{2}}(i)} = v_{i - \frac{1}{2}}},v_{i + \frac{1}{2}}$

[0038] The linear (1-th degree) interpolation on the stencil can bewritten in the Newton form as${P^{1}(v)} = {{Q\left\lbrack v_{i - \frac{1}{2}} \right\rbrack} + {{Q\left\lbrack {v_{i - \frac{1}{2}},v_{i + \frac{1}{2}}} \right\rbrack}\left( {v - v_{i - \frac{1}{2}}} \right)}}$

[0039] At the next step, since the grid is one-dimensional, only twochoices exist to expand the stencil by adding one point: either addingthe left neighbor $v_{i - \frac{3}{2}},$

[0040] resulting in the following quadratic interpolation${R(v)} = {{P^{1}(v)} = {{Q\left\lbrack {v_{i - \frac{3}{2}},v_{i - \frac{1}{2}},v_{i + \frac{1}{2}}} \right\rbrack}\left( {v - v_{i - \frac{1}{2}}} \right)\left( {v - v_{i + \frac{1}{2}}} \right)}}$

[0041] or the right neighbor $v_{i + \frac{3}{2}},$

[0042] resulting in the following quadratic interpolation${R(v)} = {{P^{1}(v)} = {{Q\left\lbrack {v_{i - \frac{1}{2}},v_{i + \frac{1}{2}},v_{i + \frac{3}{2}}} \right\rbrack}\left( {v - v_{i - \frac{1}{2}}} \right)\left( {v - v_{i + \frac{1}{2}}} \right)}}$

[0043] Note that the derivation of the above two interpolations have thesame element$\left( {v - v_{i - \frac{1}{2}}} \right)\left( {v - v_{i + \frac{1}{2}}} \right)$

[0044] multiplied by two different constants${Q\left\lbrack {v_{i - \frac{3}{2}},v_{i - \frac{1}{2}},v_{i + \frac{1}{2}}} \right\rbrack}\quad {and}\quad {{Q\left\lbrack {v_{i - \frac{1}{2}},v_{i + \frac{1}{2}},v_{i + \frac{3}{2}}} \right\rbrack}.}$

[0045] These two constants are the two second degree divided differencesof q(ν)in two different stencils. Note that a smaller divided differenceimplies the function is “smoother” in that stencil. The selection ofstencil can thus be performed by comparing the two relevant divideddifferences and picking the one with a smaller absolute value. Thus, if$\begin{matrix}{{{Q\left\lbrack {v_{i - \frac{3}{2}},v_{i - \frac{1}{2}},v_{i + \frac{1}{2}}} \right\rbrack}} < {{Q\left\lbrack {v_{i - \frac{1}{2}},v_{i + \frac{1}{2}},v_{i + \frac{3}{2}}} \right\rbrack}}} & (3)\end{matrix}$

[0046] The 3 point stencil will be selected as${{\overset{\_}{S_{3}}(i)} = v_{i - \frac{3}{2}}},v_{i - \frac{1}{2}},v_{i + \frac{1}{2}}$

[0047] Otherwise, the stencil will be selected as${{\overset{\_}{S_{3}}(i)} = v_{i - \frac{1}{2}}},v_{i - \frac{1}{2}},v_{i + \frac{3}{2}}$

[0048] The step can be continued, with one point added to the stencil ateach step, until the desired number of points in the stencil is reached.In general, to obtain a m-th order interpolating polynomial, a stencilwith m+1 points is required.

[0049] If the index grid is uniform (Δν_(i)=Δν), undivided differences$\left( {{Q\left\lbrack {v_{i + \frac{1}{2}},\cdots \quad,v_{i + j - \frac{1}{2}}} \right\rbrack} - {Q\left\lbrack {v_{i - \frac{1}{2}},\cdots \quad,v_{i + j - \frac{3}{2}}} \right\rbrack}} \right.$

[0050] can be used instead. In this case, the Newton formulae (1) and(2) for obtaining interpolating polynomials should also be adjustedaccordingly. This saves computation time and reduces round-off errors.

[0051] It is shown in [1] that for a piecewise smooth function q(ν),ENO-based method starting with a two point stencil has the followingproperties:

[0052] 1. The accuracy condition

p _(i)(ν)=q(ν)+O(Δν^(k)),ν∈I _(i)

[0053] is valid for any cell I_(i) which does not contain a C⁰ comerpoint. This implies that the ENO-based method can recover the full highorder accuracy right up to the point.

[0054] 2. p_(i)(ν) Is monotone in any cell I_(i) which does contain a“very sharp transition” or discontinuity of q(ν).

[0055] 3. The reconstruction is TVB (total variation bounded). That is,there exists a function z(ν), satisfying

z(ν)=q(ν)+O(Δν^(k−1)),ν∈I _(i)

[0056] for any cell I_(i), including those cells which containtransitions, such that

TV(z)≦TV(Q)

[0057] where TV( ) is the total variance function.

[0058] In practice, a wrong stencil may be selected due to a round-offerror of the divided difference calculation. That is, when both divideddifferences are near zero, a small change at the round off level wouldchange the direction of the inequality (1) and hence the stencilselection. To remedy this problem, a “biasing” strategy may beintroduced to weigh the selection towards a preferred stencil. In oneembodiment, the biasing strategy is implemented as follows.

[0059] First, identifying a preferred stencil${{\overset{\_}{S_{pref}}(i)} = v_{i - r - \frac{1}{2}}},\cdots \quad,v_{i - r + k + \frac{1}{2}}$

[0060] Next, replacing (3) by${{Q\left\lbrack {v_{j - \frac{1}{2}},\cdots \quad,v_{j + l - \frac{1}{2}}} \right\rbrack}} < {b{{Q\left\lbrack {v_{i + \frac{1}{2}},\cdots \quad,v_{j + l - \frac{1}{2}}} \right\rbrack}}}$if $v_{j + \frac{1}{2}} > v_{i - r + \frac{1}{2}}$

[0061] i.e. if the left-most point ν $j + \frac{1}{2}$

[0062] of the current stencil has not reached the left-most point ν$i - r + \frac{1}{2}$

[0063] of the preferred stencil; otherwise if$v_{j + \frac{i}{2}} \leq_{i - r + \frac{1}{2}}$

[0064] one replace (3) by${b{{Q\left\lbrack {v_{j - \frac{1}{2}},\cdots \quad,v_{j + l - \frac{1}{2}}} \right\rbrack}}} < {{Q\left\lbrack {v_{i + \frac{1}{2}},\cdots \quad,v_{j + l + \frac{1}{2}}} \right\rbrack}}$

[0065] Where b>1 is the biasing parameter. Analysis in [3] indicates agood choice of the parameter is b=2. The philosophy is to stay as closeas possible to the preferred stencil, unless the alternative is a factorb>1 better in smoothness.

[0066] Device models are usually described with charge and currentfunctions with three voltages as independent variables. To best reflectthe original model, the preferred embodiment uses multidimensional gridsto describe the independent variables and uses multidimensional tablesto store the grid indices and values of stencil points.

[0067] In obtaining approximation of device models in the generalm-dimensional space Ω⊂V^(m), enclosed by boundary δΩ, the most naturaland computationally efficient extension of polynomial approximations tomulti-dimension appears through the use of tensor products. FIG. 3illustrates a flow chart of this approximation method. At step 310, aset of measurement units is first identified in a (m+1)-dimensionalmanifold where m is an integer. Step 320 identifies two or more sets ofunits that neighbor one or more of said measurement units. For eachm-dimensional boundaries of the (m+1)-dimensional manifold, step 330obtains an approximation of the analytical device model based on the twosets of neighboring units and the first set of measurement units. Step340 blends the approximations obtained for m-dimensional boundaries ofthe (m+1)-dimensional manifold to produce an approximation of theanalytical device model. Steps 310-340 iterate repeatedly until the setof units achieves the desired accuracy of the analytical device model.In one embodiment, step 340 is implemented as follows.

[0068] Let D^(k), k=1, . . . ,K be K non-overlapping general hexahedralssuch that${\Omega = {\overset{K}{\bigcup\limits_{k = 1}}D^{K}}},{{\delta \quad \Omega} = {\overset{K}{\bigcup\limits_{k = 1}}{\delta \quad D_{kk}}}},{{\delta \quad D^{k}} = {\overset{K}{\bigcup\limits_{k = 1}}{\delta \quad D_{ki}}}}$

[0069] where δD_(kk)=D^(k)∩δΩ and δD_(ki)=D^(k)∩D^(i) for i≠k, and thereexists a diffeomorphism, Ψ:D^(k)→I, where I⊂V^(m) is the unit cube,i.e.I∈[−1,1]^(m). The introduction of the diffeomorphism ensures that astandard interpolation and blending scheme described in [2] can beapplied in the unit cube. The existence of the diffeomorphism isguaranteed since a circuit simulator is always free to decompose themanifold. The non-overlapping D^(k) is usually rectangular in shape withdifferent size. This is because devices have fine local behavior orsharp transitions such as the sub-threshold behavior of MOSFETS. Inorder to resolve sharp transitions and keep the table size small, anon-uniform D^(k) may be used, i.e. smaller sizes are used at sharptransition regions and bigger sizes at smooth regions.

[0070] For simplicity, D^(k) is denoted as D with boundary δD, referringto any rectangular grid table cell. The coordinates of D will be denotedas (v₁,v₂,v₂) or (ξ,η,ζ), interchangeably. DefineΞ_(L) = Ξ_(L)^(ξ) × Ξ_(L)^(η) × Ξ_(L)^(ζ)

[0071] through the subsets$\Xi_{L}^{\xi} = \left\{ {{\overset{\_}{i}\left. {\overset{\_}{i} \in \left\lbrack {0,\ldots \quad,L_{\xi}} \right\rbrack} \right\}},{\Xi_{L}^{\eta} = \left\{ {{\overset{\_}{j}\left. {\overset{\_}{j} \in \left\lbrack {0,\ldots \quad,L_{\eta}} \right\rbrack} \right\}},{\Xi_{L}^{\xi} = \left\{ {{\overset{\_}{k}\left. {\overset{\_}{k} \in \left\lbrack {0,\ldots \quad,L_{\zeta}} \right\rbrack} \right\}},} \right.}} \right.}} \right.$

[0072] Associated with these sets are nodal sets${{\Lambda_{L}^{\xi}\left( \Xi_{L}^{\xi} \right)} = \left\{ \xi_{i} \right\}_{\overset{\_}{i} = 0}^{L_{\xi}}},{{\Lambda_{L}^{\eta}\left( \Xi_{L}^{\eta} \right)} = \left\{ \eta_{\overset{\_}{j}} \right\}_{\overset{\_}{j} = 0}^{L_{\eta}}},{{\Lambda_{L}^{\zeta}\left( \Xi_{L}^{\zeta} \right)} = \left\{ \zeta_{\overset{\_}{k}} \right\}_{k = 0}^{L_{\zeta}}},$

[0073] with the global nodal set,

Λ_(L)(Ξ_(L))=Λ_(L)=Λ_(L) ^(ξ)×Λ_(L) ^(η)×Λ_(L) ^(ζ)

[0074] The nodal sets are assumed to be ordered such thatξ₀=η₀=ζ₀=−1,ξ_(L) _(ξ) =η_(L) _(η) =ζ_(L) _(ζ) =1.

[0075] Likewise, for the construction of the polynomial approximationinside the domain I, introduce the global setΞ_(N) = Ξ_(N)^(ξ) × Ξ_(N)^(η) × Ξ_(N)^(ζ) where${\Xi_{N}^{\xi} = \left\{ {{\overset{\_}{i}\left. {\overset{\_}{i} \in \left\lbrack {0,\cdots \quad,N_{\xi}} \right\rbrack} \right\}},{\Xi_{N}^{\eta} = \left\{ {\overset{\_}{j}{{\overset{\_}{j} \in \left\lbrack {0,\cdots \quad,N_{\eta}} \right\rbrack}}} \right\}},{\Xi_{N}^{\xi} = \left\{ {\overset{\_}{k}{{\overset{\_}{k} \in \left\lbrack {0,\cdots \quad,N_{\zeta}} \right\rbrack}}} \right.}} \right\}},$

[0076] with the corresponding nodal sets${{\Lambda_{N}^{\xi}\left( \Xi_{N}^{\xi} \right)} = \left\{ \xi_{i} \right\}_{\overset{\_}{i} = 0}^{N_{\xi}}},{{\Lambda_{N}^{\eta}\left( \Xi_{N}^{\eta} \right)} = \left\{ \eta_{\overset{\_}{j}} \right\}_{\overset{\_}{j} = 0}^{N_{\eta}}},{{\Lambda_{N}^{\zeta}\left( \Xi_{N}^{\zeta} \right)} = \left\{ \zeta_{\overset{\_}{k}} \right\}_{k = 0}^{N_{\zeta}}},$

[0077] For simplicity, the notationΞ_(N)^(ξ  n) = Ξ_(N)^(ξ) × Ξ_(N)^(η)

[0078] and likewise were used for various combinations of the sets.

[0079] To establish a one to one mapping between the unit cube and thegeneral rectangle-shaped cell, a global map ψ:D^(k)→I, may beconstructed in one embodiment using the so-called transfinite blendingfunctions, which are described in [2], the content of which is herebyincorporated in its entirety by reference. The map x=Ψ⁻¹ (ξ) in oneembodiment may be derived from the Boolean sumq = (P_(ξ)(q) ⊕ P_(η)(q) ⊕ P_(ζ)(q))q(ξ)   = P_(ξ)(q) + P_(η)(q) + P_(ζ)(q) − P_(ζ)P_(η)(q) − P_(ξ)P_(ζ)(q) − P_(η)P_(ζ)(q) +   P_(ξ)P_(η)P_(ζ)(q)where${P_{\xi}(q)} = {\sum\limits_{\Xi_{L}^{\xi}}{{N_{\overset{\_}{i}}^{\xi}(\xi)}{q\left( {\xi_{i},\eta,\zeta} \right)}}}$

[0080] and likewise for P_(η),(q),P_(ζ)(q), are denoted as “faceprojectors”. The “shape functions”$N_{\overset{\_}{i}}^{\xi},N_{\overset{\_}{j}}^{\eta},N_{\overset{\_}{k}}^{\zeta}$

[0081] associated with these face projectors are nothing more than theinterpolating ENO polynomials based on nodal setsΛ_(L)^(ξ), Λ_(L)^(η), Λ_(L)^(ζ)

[0082] respectively. Using the properties of projectors, the “edgeprojectors” become${P_{\xi}{P_{\eta}(q)}} = {\sum\limits_{\Xi_{L}^{\xi\eta}}{{N_{\overset{\_}{i}}^{\xi}(\xi)}{N_{\overset{\_}{j}}^{\eta}(\eta)}{q\left( {\xi_{\overset{\_}{i}},\eta_{\overset{\_}{j}},\zeta} \right)}}}$

[0083] and likewise for P_(ηζ)(q),P_(ξζ)(q). The “vertex projector”becomes${P_{\xi}P_{\eta}{P_{\zeta}(q)}} = {\sum\limits_{\Xi_{L}}{{N_{\overset{\_}{i}}^{\xi}(\xi)}{N_{\overset{\_}{j}}^{\eta}(\eta)}{N_{\overset{\_}{k}}^{\zeta}(\zeta)}{q\left( {\xi_{\overset{\_}{i}},\eta_{\overset{\_}{j}},\zeta_{\overset{\_}{k}}} \right)}}}$

[0084] To simplify things further, the following transfinite blendingfunction may also be applied to construct the face projectors.$q = {{\left( {{P_{\eta}(q)} \oplus {P_{\zeta}(q)}} \right){q\left( {\xi_{\overset{\_}{i}},\eta,\zeta} \right)}}\quad = {{P_{\eta}\left( {q\left( {\xi_{\overset{\_}{i}},\eta,\zeta} \right)} \right)} + {P_{\zeta}\left( {q\left( {\xi_{\overset{\_}{i}},\eta,\zeta} \right)} \right)} - {P_{\eta}{P_{\zeta}\left( {q\left( {\xi_{\overset{\_}{i}},\eta,\zeta} \right)} \right)}}}}$

[0085] and similarly for other face projectors. With the abovereconstruction procedure, it is guaranteed that the device modelapproximations are continuous over all the interfaces between cells.This is because they are constructed uniquely from the node points andsome neighboring points on the plane. By using the previously-mentionedbiasing scheme in ENO approximation, the derivatives are also continuousin smooth regions where the ENO approximation in each dimension of thegird chooses the points on the same side.

[0086] With the above formulations, the inverse of the global map can beconstructed in a preferred embodiment as follows. Assuming the shapefunctions are linear, i.e.

L_(ξ)=L_(η)=L_(ζ)=1,

Λ_(L) ^(ζ)=Λ_(L) ^(n)=Λ_(L) ^(ζ)={−1,1}

N ₀ ^(ξ)=(1−ξ)/2, N ₁ ^(ξ)=(1+ξ)/2

[0087] and similarly for the other shape functions:$N_{\overset{\_}{j}}^{n},N_{\overset{\_}{k}}^{\zeta}$

[0088] Thus, to construct the map a parametric form is established forthe edges enclosing D, e.g.${q\left( {\xi,{- 1},{- 1}} \right)} = {\sum\limits_{\Xi_{N}^{\xi}}{{q\left( {\xi_{i},{- 1},{- 1}} \right)}{L_{i}^{\alpha}(\xi)}}}$

[0089] Where L_(i) ^(α)(ξ) is the ENO approximating polynomial based onthe nodal set Λ_(N) _(^(ξ)) , employed in the unite cube I. For athrough discussion of transfinite blending functions and theirproperties, please refer to [2].

[0090] Though computation required in table model evaluation is just asmall percentage of analytical device model evaluation, what could causeinefficiency is the time for gathering the multi-dimensional table data.To solve this problem, one embodiment uses a cell-based data storagescheme that achieves efficient table model evaluation.

[0091] The scheme is described here assuming the grid is two-dimensional(2-D). In other embodiments, this scheme is extended to otherdimensions. The building block of 2-D table model storage schemeconsists of four neighboring cells spanning in both dimensions, as shownin FIG. 2. It may not be feasible to save all the table model entriesassociated with these four cells, since many overlapping informationwould be stored and the table size would be much bigger than necessary.As an alternative, only a part of the entries are stored to enable fastdata gathering and small table size. FIG. 4 schematically depicts aselection of the entries. The solid lines in FIG. 4 represent the ENOpolynomial information to be stored while the dashed lines representthose to be skipped. The stars represent the comer point values to bestored while the circles represent the corner points to be skipped.

[0092] With this scheme, the left-top cell contains all the informationit needs for table evaluation. The right-top and the left-bottom cellsonly need to access one neighboring building block, as shown by thearrows in FIG. 4. Only the right-bottom cell needs to access twoneighboring building blocks for information. The average is just oneneighboring building block access per cell.

[0093] In one embodiment, a dynamic table construction strategy is usedto further reduce the overlap. Specifically, when model evaluation isrequired at an operating point, only the building blocks containing theoperating point and some neighboring building blocks are constructed.This leads to big savings in memory in many cases. The building blocksmay be stored in a hash table with efficient look-up schemes for easyand fast access. Keys of the hash table may be integers representingoperating points. Values of the hash table can be build blockscontaining operating points. Hash functions can be developed usingwell-known algorithmic methods to ensure a one-to-one mapping betweenthe keys and values of the hash table and efficient look-up.

[0094] The building blocks or grid cells may be of either uniform ornon-uniform size. For evaluation of electronic device models, the gridis non-uniform in nature since electronic devices, such as transistors,may exhibit sharp transition behavior in certain grid regions and finercells are required to reflect these sharp transitions. FIG. 5schematically depicts a non-uniform grid for approximating the behaviorof a MOSFET transistor. The solid lines represent the edges that needapproximation. For some edges of the smaller-sized cells, there is nogrid point to the left or the right in the table, the dashed lines showthe extended, but not stored, points that is needed to build thenon-uniform table. In one embodiment, the circles are considered as partof the stencil in performing ENO approximation.

[0095] In one embodiment, it possible to combine ENO approximating ofdifferent orders is combined on different edges. For example, athird-order approximation in the exponential region of BJTs while usingquadratic or even linear approximation in other regions may be applied.With the use of transfinite blending functions, the approximation isstill smooth across shared boundaries.

[0096] One embodiment has been tested in the simulation of thousands ofdigital, analog, and mixed-signal circuits. Results show that the methoddescribed above turns out to be very robust and does not haveconvergence problems. FIG. 6 shows the speed improvement of using theENO-based table model method described above in the simulation of aswitched-capacitor circuit. By using the ENO-based method, a speed-up of5 times over the same simulation using the prior art analytical modelsimulation has been achieved. The actual simulation data is shown inTable 1. associated with FIG. 6. This fact shows that the ENO-basedmethod is also more accurate than the prior art analytical modelsimulation.

[0097] The speed and accuracy of the ENO-based table model methoddescribed above can be further demonstrated by additional simulations.FIG. 7 shows the result of the simulation of a sigma-delta convertercircuit containing 12222 MOSEFTs and hundreds of other devices. TheENO-based method yields a speed-up of 3.5 times, as shown by Table 1.One can see from FIG. 7 that there is almost no difference between theapproximation and the actual curve. FIG. 8 demonstrates the simulationof a digital-analog converter containing 12553 MOSEFTs using theENO-based method. The result shows that the ENO-based method matches theoriginal curve even for the initial oscillations. The simulation is 4.5times faster than prior art methods that use analytical model. FIG. 9shows the simulation of a memory circuit consisting of 108267 elements,most of which are MOSEFTs. The ENO-based method is more than two timesfaster than prior art methods that use analytical model. Finally, FIG.10 shows simulation of a clock-tree circuit consisting of 2833737elements, most of which are MOSEFTs, diodes, and RC parasitic. TheENO-based method finishes the transient simulation and delaymeasurements in just a couple of thousands of seconds. TABLE 1Simulation Results in Seconds Model-Test Switch-Cap Filter Sigma-DeltaD/A Converter ENO table model 1580 8900 7150 Analytical model 8072 3300032100

[0098] While the above invention has been described with reference tocertain preferred embodiments, the scope of the present invention is notlimited to these embodiments. One skilled in the art may find variationsof these preferred embodiments which, nevertheless, fall within thespirit of the present invention, whose scope is defined by the claimsset forth below.

[0099] References

[0100] [1] A. Harten, B. Engquist, S. Osher and S. Chakravarthy, “Someresults on uniformly high order accurate essentially non-oscillatoryschemes,” Applied Numerical Mathematics, v2 (1986), pp. 347-377.

[0101] [2] W. J. Gorgon and C. A. Hall, “Transfinite element methods:Blending-Function interpolation over arbitrary curved element domains,”Numerical Mathematics, vol. 21, pp. 109-129, 1973.

[0102] [3] Chi-Wang Su, “Numerical experiments on the accuracy of ENOand modified ENO schemes,” Journal of Scientific Computing, v5 (1990),pp. 127-149.

What is claimed is:
 1. A method for obtaining smooth and accurateapproximations of analytical device models, comprising the steps of:identifying a first set of measurement units; locating two or more setsof units that neighbor one or more of said measurement units; for eachof the two or more sets of neighbor units, obtaining the union of one ormore of said sets of neighbor units and said first set of measurementunits; calculating the smoothness of the analytical device model withinone or more of said unions; and selecting the union within which theanalytical device model is the smoothest as a new set of measurementunits.
 2. The method of claim 1, wherein the step of selecting furthercomprises the steps of: selecting a preferred one of said unions; andweighting the smoothness towards said preferred union.
 3. The method ofclaim 1, further comprising the step of: recording the calculatedsmoothness of the analytical device model within the selected union. 4.The method of claim 1, further comprising the step of: constructing anapproximation of the analytical device model based on the new set ofmeasurement units and the smoothness calculation.
 5. The method of claim4, wherein the approximation is a polynomial approximation.
 6. Themethod of claim 1, wherein the first set of measurement units are a setof adjacent cells on a grid.
 7. The method of claim 6, wherein theneighboring units are located based on cells neighboring the set ofadjacent cells on the grid.
 8. The method of claim 6, further comprisingthe step of: reducing the size of the cells on the grid to produce arefined set of stencil points.
 9. The method of claim 4, wherein thefirst set of measurement units are a set of adjacent cells${I_{i} \equiv \left\lbrack {v_{i - \frac{1}{2}},v_{i + \frac{1}{2}}} \right\rbrack},{i = 1},2,{\ldots \quad N}$

on a one-dimensional grid$a = {v_{\frac{1}{2}} < v_{\frac{3}{2}} < \quad {\ldots \quad v_{N - \frac{1}{2}}} < {v_{N + \frac{1}{2}}.}}$


10. The method of claim 9, wherein the neighboring units are cells$\left\lbrack {v_{i - \frac{3}{2}},v_{i - \frac{1}{2}}} \right\rbrack \quad {{{and}\quad\left\lbrack {v_{i + \frac{1}{2}},v_{i + \frac{3}{2}}} \right\rbrack}.}$


11. The method of claim 10, wherein the smoothness of the analyticaldevice model within each union of units is calculated by:${{Q\left\lbrack v_{i - r - \frac{1}{2}} \right\rbrack} \equiv {q\left( v_{i - r - \frac{1}{2}} \right)}},{i = 1},2,{\cdots \quad {N.{and}}}$${Q\left\lbrack {v_{i - r - \frac{1}{2}},\cdots \quad,v_{i - r + j - \frac{1}{2}}} \right\rbrack} \equiv {\frac{{Q\left\lbrack {v_{i - r + \frac{1}{2}},\cdots \quad,v_{i - r + j - \frac{1}{2}}} \right\rbrack} - {Q\left\lbrack {v_{i - r - \frac{1}{2}},\cdots \quad,v_{i - r + j - \frac{3}{2}}} \right\rbrack}}{v_{i - r + j - \frac{1}{2}} - v_{i - r - \frac{1}{2}}}{for}\quad j} \geq 1$

where r is either 0 or
 1. 12. The method of claim 10, wherein theapproximation of the analytical model is constructed by:${p(v)} = {\sum\limits_{j = 0}^{k}\quad {{Q\left\lbrack {v_{i - r - \frac{1}{2}},\cdots \quad,v_{i - r + j - \frac{1}{2}}} \right\rbrack}{\sum\limits_{{m = 0},{m \neq r}}^{j - 1}\quad {\prod\limits_{m = 0}^{j - 1}\quad \left( {v - v_{i - r + m - \frac{1}{2}}} \right)}}}}$

where r is either 0 or
 1. 13. A method for obtaining smooth and accurateapproximations of analytical device models in a m-dimensional manifold,comprising the steps of: providing a partition of the measurement space,and for each partition, executing the steps of claim
 1. 14. A method forobtaining smooth and accurate approximations of analytical device modelsin a (m+1)-dimensional manifold where m is an integer, comprising thesteps of: identifying a first set of measurement units in the(m+1)-dimensional manifold; locating two or more set of units thatneighbor one or more of said measurement units; obtaining, for eachm-dimensional boundary of the (m+1)-dimensional manifold, anapproximation of the analytical device model based on the two or moresets of neighboring units and the first set of measurement units; andblending the approximations obtained for each m-dimensional boundary ofthe (m+1)-dimensional manifold to produce an approximation of theanalytical device model.
 15. The method of claim 14, wherein the step oflocating further comprises the step of: locating two or more sets ofunits in a subspace of the (m+1)-dimensional manifold.
 16. The method ofclaim 14, wherein the step of obtaining further comprises the steps of:obtaining, for each m-dimensional boundary of the (m+1)-dimensionalmanifold, a union of each set of the two or more sets of neighboringunits and said first set of measurement units; calculating thesmoothness of the approximation for each said union; and selecting theunion that renders the smoothest approximation of the analytical devicemodel.
 17. The method of claim 16, wherein the step of obtaining furthercomprises the step of: mapping each of the two or more sets ofneighboring units and the first set of measurement units to a standardunit element in the (m+1)-dimensional manifold.
 18. The method of claim16, wherein the (m+1)-dimensional manifold is a (m+1)-dimensional grid.19. The method of claim 16, wherein the standard unit element is a unitcube.
 20. The method of claim 17, wherein the step of blending furthercomprises the steps of: projecting the faces, edges, and vertices ofeach standard unit element to an approximated curve based on theapproximations obtained for each m-dimensional boundary of the(m+1)-dimensional manifold.
 21. The method of claim 14, wherein the stepof identifying further comprises the step of identifying a first set ofmeasurement units in the neighborhood of an operating point of theanalytical device model.
 22. A method for reducing memory usage duringcircuit simulation, comprising the steps of: identifying one or moreoperating points; obtaining, for each of the one or more operatingpoints, a set of measurement points surrounding said operating point;eliminating overlapping boundaries among each set of measurement points;and constructing a hash table for storing each set of measurementpoints.
 23. The method of claim 22, wherein said one or more operatingpoints are identified dynamically during circuit simulation.
 24. Themethod of claim 22, wherein the step of constructing a hash tablefurther comprises the steps of: transforming each of the one or moreoperating points into a unique integer; and constructing a hashfunction, said hash function providing a one-to-one mapping between theunique integers and the sets of measurement points.